The EEAaq package allows users to retrieve air quality data for multiple geographical zones, pollutants, and time periods in a single request. Queries are submitted as lists, which enables flexibility in specifying combinations of parameters.

library(EEAaq)
library(tidyverse)

## Warning: il pacchetto 'ggplot2' è stato creato con R versione 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

EEAaq_get_data

Below we demonstrate the use of query by using different combinations of user-defined arguments.

Retrieve NO\(_2\) data for a specific municipality (LAU zone) given its unique identifier `LAU_ID`

data_lau <- EEAaq::EEAaq_get_data(
  zone_name = "15146",      # LAU zone code
  NUTS_level = "LAU",       # NUTS level
  LAU_ISO = "IT",           # Country code for Italy
  pollutants = "PM10",      # Pollutant 
  from = "2022-01-01",      # Start date
  to = "2023-12-31",        # End date
  verbose = FALSE           # Print detailed progress
)

# Preview the first few rows of the dataset
head(data_lau)

Retrieve NO\(_2\) data for a specific macroregion (Eurostat classification NUTS-1) given its name `LATN_NAME`

# Identify the names of the areas from which to download the data
zones <- c("Région de Bruxelles-Capitale/Brussels Hoofdstedelijk Gewest","Vlaams Gewest","West-Nederland","Zuid-Nederland")

# Download the corresponding data
data <- EEAaq_get_data(
  zone_name = zones,                 # LAU zone code
  NUTS_level = "NUTS1",              # NUTS level
  pollutants = c("NO2", "PM10"),     # Pollutant 
  from = "2023-01-01",               # Start date
  to = "2023-12-31",                 # End date
  verbose = FALSE                    # Print detailed progress
)

unique(data$AirQualityStationEoICode)

##  [1] "BELAL01" "BELAT83" "BELHB23" "BETB001" "BETB004" "BETB006" "BETB008"
##  [8] "BETB011" "BETBUL1" "BETCHA1" "BETE013" "BETE714" "BETE716" "BETM802"
## [15] "BETMEU1" "BETN043" "BETR001" "BETR002" "BETR012" "BETR701" "BETR702"
## [22] "BETR721" "BETR740" "BETR801" "BETR802" "BETR803" "BETR804" "BETR805"
## [29] "BETR806" "BETR817" "BETR818" "BETR822" "BETR831" "BETR842" "BETR891"
## [36] "BETR897" "BETREG1" "BETVBX1" "BETVBX2" "BETVBX3" "NL00136" "NL00138"
## [43] "NL00236" "NL00237" "NL00240" "NL00241" "NL00247" "NL00546" "NL00551"
## [50] "NL00553" "NL00556" "NL00570" "NL00572" "NL00573" "NL00701" "NL00704"

Note 1: If the query’s zone_name parameter corresponds to a valid CITY_NAME (i.e., not NULL in the dataset), the function will return the corresponding data. If no valid CITY_NAME is associated with the zone_name, the function attempts to retrieve all available data for the entire country and subsequently filter for the specified zone_name.

Note 2: For very small towns or certain countries such as Turkey or Albania, data may not currently be available in the dataset.This limitation reflects the data unavailability at the EEA Air Quality Viewer.

Note 3: If the parameters used in the query include polygon or quadrant, the function outputs an EEAaq_df_sfc object. Otherwise, it returns an EEAaq_df object, which is a tibble dataframe.

EEAaq map stations

EEAaq_map_stations generates a static or dynamic map of user-defined monitoring stations. The function accepts as input either an object of the EEAaq_df class (default output of the EEAaq_get_data function), or all other parameters specifying the area and the pollutants.

Map the stations using as `EEAaq_df` object, the dataset concerning NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands

EEAaq_map_stations(
  data = data,
  bounds_level = "NUTS3",
  color = FALSE,
  dynamic = FALSE
)

## Simple feature collection with 4 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2.546088 ymin: 50.688 xmax: 6.225231 ymax: 53.18511
## Geodetic CRS:  WGS 84
##   NUTS_ID LEVL_CODE CNTR_CODE
## 1     BE1         1        BE
## 2     BE2         1        BE
## 3     NL3         1        NL
## 4     NL4         1        NL
##                                                     NAME_LATN
## 1 Région de Bruxelles-Capitale/Brussels Hoofdstedelijk Gewest
## 2                                               Vlaams Gewest
## 3                                              West-Nederland
## 4                                              Zuid-Nederland
##                                                     NUTS_NAME MOUNT_TYPE
## 1 Région de Bruxelles-Capitale/Brussels Hoofdstedelijk Gewest         NA
## 2                                               Vlaams Gewest         NA
## 3                                              West-Nederland         NA
## 4                                              Zuid-Nederland         NA
##   URBN_TYPE COAST_TYPE                       geometry
## 1        NA         NA MULTIPOLYGON (((4.415738 50...
## 2        NA         NA MULTIPOLYGON (((5.776583 50...
## 3        NA         NA MULTIPOLYGON (((5.171192 52...
## 4        NA         NA MULTIPOLYGON (((5.518671 51...
## points Country ISO AirQualityStationEoICode AirQualityStationNatCode AirQualityStationName Altitude NUTS1 NUTS1_ID NUTS2 NUTS2_ID NUTS3 NUTS3_ID LAU_NAME LAU_ID AirPollutant geometry n

Note: Using the parameter bounds_level = "NUTS3", the map is generated with internal boundaries corresponding to the NUTS-3 level. The same output could be obtained specifying explicitly the zone information.

Map all the stations monitoring NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands

EEAaq_map_stations(
  zone_name = zones,
  NUTS_level = "NUTS1",
  pollutant = c("NO2", "PM10"),
  bounds_level = "NUTS3",
  color = FALSE,
  dynamic = FALSE
)

## Simple feature collection with 4 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2.546088 ymin: 50.688 xmax: 6.225231 ymax: 53.18511
## Geodetic CRS:  WGS 84
##   NUTS_ID LEVL_CODE CNTR_CODE
## 1     BE1         1        BE
## 2     BE2         1        BE
## 3     NL3         1        NL
## 4     NL4         1        NL
##                                                     NAME_LATN
## 1 Région de Bruxelles-Capitale/Brussels Hoofdstedelijk Gewest
## 2                                               Vlaams Gewest
## 3                                              West-Nederland
## 4                                              Zuid-Nederland
##                                                     NUTS_NAME MOUNT_TYPE
## 1 Région de Bruxelles-Capitale/Brussels Hoofdstedelijk Gewest         NA
## 2                                               Vlaams Gewest         NA
## 3                                              West-Nederland         NA
## 4                                              Zuid-Nederland         NA
##   URBN_TYPE COAST_TYPE                       geometry
## 1        NA         NA MULTIPOLYGON (((4.415738 50...
## 2        NA         NA MULTIPOLYGON (((5.776583 50...
## 3        NA         NA MULTIPOLYGON (((5.171192 52...
## 4        NA         NA MULTIPOLYGON (((5.518671 51...
## points Country ISO AirQualityStationEoICode AirQualityStationNatCode AirQualityStationName Altitude NUTS1 NUTS1_ID NUTS2 NUTS2_ID NUTS3 NUTS3_ID LAU_NAME LAU_ID AirPollutant geometry n

Dynamic (interactive) map of all the PM\(_{10}\) monitoring stations in Milano

EEAaq_map_stations(
  data = data_lau,
  color = TRUE,
  dynamic = TRUE,
  ID = TRUE
)

## Simple feature collection with 1 feature and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 9.042173 ymin: 45.38829 xmax: 9.277104 ymax: 45.53535
## Geodetic CRS:  WGS 84
##    GISCO_ID ISO LAU_ID LAU_NAME POP_2021 POP_DENS_2 AREA_KM2 YEAR NUTS3_ID
## 1 IT_015146  IT  15146   Milano  1397715         NA 181.8155 2021    ITC4C
##                         geometry
## 1 MULTIPOLYGON (((9.120405 45...
## points Country ISO AirQualityStationEoICode AirQualityStationNatCode AirQualityStationName Altitude NUTS1 NUTS1_ID NUTS2 NUTS2_ID NUTS3 NUTS3_ID LAU_NAME LAU_ID AirPollutant geometry n

EEAaq summary

This function aims to describe the dataset that has been previously imported, both at a global level, which means considering the complete set of time stamps and monitoring stations in the dataset, and at the station-specific level, where summary statistics and information are grouped by monitoring station.
In addition to basic exploratory descriptive statistics (e.g., average pollutant concentration, variability, measures of skewness and kurtosis), the function provides information about the gap length and the correlation between pollutants if at least two pollutants are considered simultaneously.
The EEAaq_summary function receives as input an EEAaq_df object, i.e. the output of the EEAaq get data function.

Compute the descriptive statistics

summ <- EEAaq_summary(data = data)

## The dataset contains:
##  ** 477510 total observations 
##  ** 56 stations 
##  ** 8736 time stamps: from 2023-01-01 01:00:00 to 2023-12-31

Print screen the global statsitics

summ$Summary

Print screen the station-specific statsitics

summ$Summary_byStat$Mean_byStat

Print screen the linear correlation matrix

summ$Corr_Matrix

EEAaq time aggregate

Recall that most pollutants are monitored by EEA on a hourly or daily basis, posing challenges for interpretation and representation. The EEAaq_time_aggregate function simplifies this by aggregating data into annual, monthly, weekly, daily, or hourly intervals, generating summary statistics for each station in an EEAaq_taggr_df object.

Get the station-specific monthly minimum, maximum, average and median concentrations of NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands

t_aggr <- EEAaq_time_aggregate(
  data = data,
  frequency = "monthly",
  aggr_fun = c("min", "max", "mean", "median" )
)

Print screen of the aggregated (monthly) data

t_aggr$TimeAggr

Print screen of the PM\(_{10}\) aggregated data only

t_aggr$TimeAggr_byPollutant$PM10

EEAaq_idw_map

To enable quick and intuitive visual analysis, the EEAaq_idw_map function provides spatial interpolation maps using the Inverse Distance Weighting (IDW) method (Shepard, 1968). This technique estimates the value of a variable at unknown locations by calculating a weighted average of known values, with weights inversely proportional to the distance from known points. Closer points contribute more heavily to the estimate, making it a practical approach for interpolating geolocated air quality data.

Generate IDW interpolated maps of monthly average concentrations of NO\(_2\) in the Netherlands and Belgium

EEAaq::EEAaq_idw_map(
  data = t_aggr,
  pollutant = "PM10",
  aggr_fun = "mean",
  distinct = TRUE,
  gradient = TRUE,
  idp = 2
)

## Map initialization started at 2025-02-05 13:30:08.525219
## Map initialization ended at 2025-02-05 13:30:19.841176
## Computing IDW interpolation started at 2025-02-05 13:30:19.841352
## Computing IDW interpolation for: 2023-01-01, 1 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-02-01, 2 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-03-01, 3 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-04-01, 4 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-05-01, 5 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-06-01, 6 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-07-01, 7 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-08-01, 8 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-09-01, 9 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-10-01, 10 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-11-01, 11 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-12-01, 12 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation ended at 2025-02-05 13:31:54.439807

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

Generate IDW interpolated maps of the maximum monthly concentrations of NO\(_2\) in january and february 2023 in the Netherlands and Belgium

EEAaq::EEAaq_idw_map(
  data = t_aggr$TimeAggr_byPollutant$PM10 %>% dplyr::filter(Date %in% c("2023-01-01","2023-02-01")),
  pollutant = "PM10",
  aggr_fun = "max",
  distinct = TRUE,
  gradient = TRUE,
  idp = 2
)

## Map initialization started at 2025-02-05 13:32:01.657422
## Map initialization ended at 2025-02-05 13:32:14.605266
## Computing IDW interpolation started at 2025-02-05 13:32:14.605443
## Computing IDW interpolation for: 2023-01-01, 1 of 2
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-02-01, 2 of 2
## [inverse distance weighted interpolation]
## Computing IDW interpolation ended at 2025-02-05 13:32:29.484387

## [[1]]

## 
## [[2]]

Generate IDW interpolated maps of monthly average concentrations of PM\(_{10}\) in the Netherlands and Belgium

EEAaq::EEAaq_idw_map(
  data = t_aggr,
  pollutant = "PM10",
  aggr_fun = "mean",
  distinct = TRUE,
  gradient = FALSE,
  fill_NUTS_level = "NUTS3",
  dynamic = TRUE
)

## Map initialization started at 2025-02-05 13:32:30.684601
## Map initialization ended at 2025-02-05 13:32:42.167963
## Computing IDW interpolation started at 2025-02-05 13:32:42.16821
## Computing IDW interpolation for: 2023-01-01, 1 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-02-01, 2 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-03-01, 3 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-04-01, 4 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-05-01, 5 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-06-01, 6 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-07-01, 7 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-08-01, 8 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-09-01, 9 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-10-01, 10 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-11-01, 11 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-12-01, 12 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation ended at 2025-02-05 13:34:23.277381

Note: By setting fill_NUTS_level = "NUTS3", we require the predictions/interpolations to be aggregated (averaged) at the NUTS-3 level.

Extended vignette for EEAaq: Handle Air Quality Data from the European Environment Agency Data Portal

Extended vignette

Paolo Maranzano, University of Milano-Bicocca, Italy, paolo.maranzano@unimib.it

Riccardo Borgoni, University of Milano-Bicocca, Italy, riccardo.borgoni@unimib.it

Agostino Tassan Mazzocco, University of Milano-Bicocca, Italy

Samir Doghmi, University of Milano-Bicocca, Italy

05 febbraio 2025

EEAaq_get_data

Retrieve NO\(_2\) data for a specific municipality (LAU zone) given its unique identifier `LAU_ID`

Retrieve NO\(_2\) data for a specific macroregion (Eurostat classification NUTS-1) given its name `LATN_NAME`

EEAaq map stations

Map the stations using as `EEAaq_df` object, the dataset concerning NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands

Map all the stations monitoring NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands

Dynamic (interactive) map of all the PM\(_{10}\) monitoring stations in Milano

EEAaq summary

Compute the descriptive statistics

Print screen the global statsitics

Print screen the station-specific statsitics

Print screen the linear correlation matrix

EEAaq time aggregate

Get the station-specific monthly minimum, maximum, average and median concentrations of NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands

Print screen of the aggregated (monthly) data

Print screen of the PM\(_{10}\) aggregated data only

EEAaq_idw_map

Generate IDW interpolated maps of monthly average concentrations of NO\(_2\) in the Netherlands and Belgium

Generate IDW interpolated maps of the maximum monthly concentrations of NO\(_2\) in january and february 2023 in the Netherlands and Belgium

Generate IDW interpolated maps of monthly average concentrations of PM\(_{10}\) in the Netherlands and Belgium

Extended vignette for EEAaq: Handle Air Quality Data from the European Environment Agency Data Portal

Extended vignette

Paolo Maranzano, University of Milano-Bicocca, Italy, paolo.maranzano@unimib.it

Riccardo Borgoni, University of Milano-Bicocca, Italy, riccardo.borgoni@unimib.it

Agostino Tassan Mazzocco, University of Milano-Bicocca, Italy

Samir Doghmi, University of Milano-Bicocca, Italy

05 febbraio 2025

EEAaq_get_data

Retrieve NO\(_2\) data for a specific municipality (LAU zone) given its unique identifier LAU_ID

Retrieve NO\(_2\) data for a specific macroregion (Eurostat classification NUTS-1) given its name LATN_NAME

EEAaq map stations

Map the stations using as EEAaq_df object, the dataset concerning NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands

Map all the stations monitoring NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands

Dynamic (interactive) map of all the PM\(_{10}\) monitoring stations in Milano

EEAaq summary

Compute the descriptive statistics

Print screen the global statsitics

Print screen the station-specific statsitics

Print screen the linear correlation matrix

EEAaq time aggregate

Get the station-specific monthly minimum, maximum, average and median concentrations of NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands

Print screen of the aggregated (monthly) data

Print screen of the PM\(_{10}\) aggregated data only

EEAaq_idw_map

Generate IDW interpolated maps of monthly average concentrations of NO\(_2\) in the Netherlands and Belgium

Generate IDW interpolated maps of the maximum monthly concentrations of NO\(_2\) in january and february 2023 in the Netherlands and Belgium

Generate IDW interpolated maps of monthly average concentrations of PM\(_{10}\) in the Netherlands and Belgium

Retrieve NO\(_2\) data for a specific municipality (LAU zone) given its unique identifier `LAU_ID`

Retrieve NO\(_2\) data for a specific macroregion (Eurostat classification NUTS-1) given its name `LATN_NAME`

Map the stations using as `EEAaq_df` object, the dataset concerning NO\(_2\) and PM\(_{10}\) in Belgium and The Netherlands